Three-dimensional elastic frequency-domain iterative solver for full waveform inversion

ABSTRACT

Certain implementations of a three-dimensional elastic frequency domain iterative solver for full waveform inversion can be implemented as a method in which frequency domain numerical simulation of elastic waves is per formed in three-dimensional (3D) media.

TECHNICAL FIELD

This disclosure relates to methods for numerical simulation of elasticwavefields. More particularly, the disclosure relates to methods forsolving three dimensional isotropic elastic wave equation in thetemporal frequency domain for application in Full-waveform inversion(FWI).

BACKGROUND

Full-waveform inversion (FWI) is a tool for deriving velocity models inareas with complex geological conditions. FWI may be implemented in thetime domain or in the frequency domain (FD FWI). Each realization hassome advantages and disadvantages. For example, FD FWI can allowobtaining satisfactory results by hierarchical inversion of a smallgroup of low frequencies. One difficulty of FD FWI is the lack of anefficient frequency domain solver. Ideally, the solver should beoriented for situations which may be described briefly as:3D+elasticity+anisotropy+attenuation. Computing wave, even in isotropiccases, fields using a 3D frequency-domain solver in such a case can becomplicated because huge linear systems must be solved. In suchinstances, iterative solvers can be more beneficial compared to directsolvers, for example, due to lesser memory requirements.

SUMMARY

This disclosure describes technologies relating to numerical simulationof elastic wavefields, for example, using an iterative method. Thetechniques described in this disclosure can be implemented as athree-dimensional (3D) elastic frequency-domain iterative solver forFWI.

Certain implementations of the subject matter described here can beimplemented as a method including performing frequency-domain numericalsimulation of elastic waves in three-dimensional (3D) media. The 3Dmedia can include heterogeneous media. The 3D media can includeisotropic media. The 3D media can include a rock formation. With orwithout any one or more of the preceding features, the 3D media caninclude a hydrocarbon reservoir. With or without any one or more of thepreceding features, performing the frequency-domain numerical simulationcan include representing the 3D elastic system in a compliance form.With or without any one or more of the preceding features, representingthe 3D elastic system in a compliance form can include representing the3D elastic system by the following equation:

${Lu} = {{\left\lbrack {{i\omega {M\left( {x,y,z} \right)}} - {P{\frac{\partial}{\partial x}{- Q}}\frac{\partial}{\partial y}} - {R\frac{\partial}{\partial z}}} \right\rbrack u} = {f.}}$

With or without any one or more of the preceding features, representingthe 3D elastic system in a compliance form can include settingparameters of the 3D media, dimensions of a computational domain and acondition of the top boundary. With or without any one or more of thepreceding features, setting parameters of the 3D media, dimensions of acomputational domain and a condition of the top boundary comprisesrepresenting to parameters included in the representation of the 3Delastic system by the following equations:

${M = \begin{pmatrix}{\rho \; I_{3 \times 3}} & 0 \\0 & S_{6 \times 6}\end{pmatrix}},{S = \begin{bmatrix}a & {- b} & {- b} & 0 & 0 & 0 \\{- b} & a & {- b} & 0 & 0 & 0 \\{- b} & {- b} & a & 0 & 0 & 0 \\0 & 0 & 0 & c & 0 & 0 \\0 & 0 & 0 & 0 & c & 0 \\0 & 0 & 0 & 0 & 0 & c\end{bmatrix}},{a = \frac{\lambda + \mu}{\mu \left( {{2\mu} + {3\lambda}} \right)}},{b = \left( \frac{\lambda}{2{\mu \left( {{2\mu} + {3\lambda}} \right)}} \right)},{c = {\left( \frac{1}{\mu} \right).}}$

With or without any one or more of the preceding features, performingthe frequency-domain numerical simulation can include calculatingone-dimensional (1D) background medium parameters and constructing thepreconditioner. With or without any one or more of the precedingfeatures, constructing the preconditioner can include constructing thepreconditioner as an inverse to the operator shown in the followingequation:

$L_{0} = {{i\omega {M_{0}(z)}} - {P\frac{\partial}{\partial x}} - {Q\frac{\partial}{\partial y}} - {R{\frac{\partial}{\partial z}.}}}$

With or without any one or more of the preceding features, the operatorM₀ has the same structure as the operator M. With or without any one ormore of the preceding features, a complex damping function is introducedinto the operator M₀. With or without any one or more of the precedingfeatures, the complex damping function can be represented by thefollowing equations: a₀(z)={circumflex over (β)}â₀(z), b₀(z)={circumflexover (β)}{circumflex over (b)}₀(z), c₀(z)={circumflex over (β)}ĉ₀(z).{circumflex over (β)}=(1+iβ), 0<β<1. With or without any one or more ofthe preceding features, performing the frequency-domain numericalsimulation can include representing the initial problem as theperturbation to the preconditioner. With or without any one or more ofthe preceding features, representing the initial problem as theperturbation to the preconditioner can include representing the initialoperator L as a combination of L₀ and δL. With or without any one ormore of the preceding features, representing the initial operator L as acombination of L₀ and δL can include representing the initial operatoraccording to the following equation: δL=−iω[M₀(z)−M(x,y,z)]. With orwithout any one or more of the preceding features, performing thefrequency-domain numerical simulation can include decomposing thecomputational domain in one of the lateral direction. With or withoutany one or more of the preceding features, performing thefrequency-domain numerical simulation comprises starting an iterativeprocess. With or without any one or more of the preceding features, theiterative process can include a parallel version of biconjugate gradientstabilized (BiCGSTAB) method. With or without any one or more of thepreceding features, performing the frequency-domain numerical simulationcan include performing matrix-vector multiplication at each iteration ofBiCGSTAB. With or without any one or more of the preceding features, thematrix-vector multiplication can be performed by the application of aparallel Fast Fourier Transform (FFT), a banded matrix solver, and aninverse FFT. With or without any one or more of the preceding features,performing the frequency-domain numerical simulation comprisesmultiplying the solution by an inverse operator to the preconditioner.

Certain aspects of the subject matter described here can be implementedas a computer-readable medium storing instructions executable by dataprocessing apparatus to perform operations described here.

Certain aspects of the subject matter described here can be implementedas a system that includes a data processing apparatus, and acomputer-readable medium storing instructions executable by the dataprocessing apparatus to perform operations described here.

The details of one or more implementations of the subject matterdescribed in this specification are set forth in the accompanyingdrawings and the description below. Other features, aspects, andadvantages of the subject matter will become apparent from thedescription, the drawings, and the claims.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flowchart of an example of a process for numericalsimulation of 3D elastic waves in frequency domain.

FIG. 2 is a schematic diagram representing message passing interface(MPI) parallelization with domain decomposition in a Y-direction.

FIG. 3 is a plot comparing simulated horizontal velocity and atime-domain solver for a 3D heterogeneous model with absorbing boundarycondition at the top.

FIG. 4 is a plot comparing simulated horizontal velocity and atime-domain solver for a 3D heterogeneous model with free surfaceboundary condition at the top.

FIG. 5 is a plot showing scalability of a parallel iterative solver.

FIG. 6 is a plot showing velocity distribution in a land model.

FIG. 7 is a schematic diagram showing a real part of a frequency-domainhorizontal velocity component.

DETAILED DESCRIPTION

This disclosure describes an iterative method for frequency-domainnumerical simulation of elastic waves in three-dimensional (3D)heterogeneous, isotropic, land media. The convergence of the numericalsimulation is achieved by implementing a preconditioner tuned for 3Dsimulations. The techniques described here can be implemented for fullwaveform inversion (FWI). For example, the techniques described here canbe implemented at low frequencies and in regions with more verticalvariations than lateral variations. Some implementations of the subjectmatter described here can be implemented using computing clusters, forexample, high performance computing clusters with distributed memory.The algorithm described here shows good and strong scalability. Suchimplementations allow efficient simulation in target 3D media of largesize.

In some implementations, the techniques described in this disclosure canbe used to simulate 3D elastic waves in frequency domain.

As a first step, the 3D elastic system is written in compliance form asshown in Equation (1) below.

$\begin{matrix}{{Lu} = {{\left\lbrack {{i\; \omega \; {M\left( {x,y,z} \right)}} - {P\frac{\partial}{\partial x}} - {Q\frac{\partial}{\partial y}} - {R\frac{\partial}{\partial z}}} \right\rbrack u} = f}} & (1)\end{matrix}$

In Equation (1), M is an operator represented by Equation (2) below.

$\begin{matrix}{M = \begin{pmatrix}{\rho \; I_{3 \times 3}} & 0 \\0 & S_{6 \times 6}\end{pmatrix}} & (2)\end{matrix}$

In Equation (2), ρ is density, I_(3×3) is an identical operator andS_(6×6) is a matrix shown in Equation (3) below.

$\begin{matrix}{S = \begin{bmatrix}a & {- b} & {- b} & 0 & 0 & 0 \\{- b} & a & {- b} & 0 & 0 & 0 \\{- b} & {- b} & a & 0 & 0 & 0 \\0 & 0 & 0 & c & 0 & 0 \\0 & 0 & 0 & 0 & c & 0 \\0 & 0 & 0 & 0 & 0 & c\end{bmatrix}} & (3)\end{matrix}$

Also, in Equation (2), a, b, and c are operators represented byEquations (4)-(6) below.

$\begin{matrix}{a = \frac{\lambda + \mu}{\mu \left( {{2\mu} + {3\lambda}} \right)}} & (4) \\{b = \left( \frac{\lambda}{2{\mu \left( {{2\mu} + {3\lambda}} \right)}} \right)} & (5) \\{c = \left( \frac{1}{\mu} \right)} & (6)\end{matrix}$

In Equations (4)-(6), λ and μ are Lame parameters. In addition,S_(6×6)=C⁻¹, a compliance tensor.

Returning to Equation (1), P, Q and R are represented by Equations(7)-(9) below.

$\begin{matrix}{P = \begin{bmatrix}0 & A \\A^{T} & 0\end{bmatrix}} & (7) \\{Q = \begin{bmatrix}0 & B \\B^{T} & 0\end{bmatrix}} & (8) \\{R = \begin{bmatrix}0 & D \\D^{T} & 0\end{bmatrix}} & (9)\end{matrix}$

A, B and D in Equations (7)-(9) are represented by Equations (10)-(12).

$\begin{matrix}{A = \begin{bmatrix}1 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 1 \\0 & 0 & 0 & 0 & 1 & 0\end{bmatrix}} & (10) \\{B = \begin{bmatrix}0 & 0 & 0 & 0 & 0 & 1 \\0 & 1 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 1 & 0 & 0\end{bmatrix}} & (11) \\{D = \begin{bmatrix}0 & 0 & 0 & 0 & 1 & 0 \\0 & 0 & 0 & 1 & 0 & 0 \\0 & 0 & 1 & 0 & 0 & 0\end{bmatrix}} & (12)\end{matrix}$

Also, in Equation (1), u=(u_(x), u_(y), u_(z), σ_(xx), σ_(yy), σ_(zz),σ_(xz), σ_(xy)), f(x _(s), ω) determines the type of source and itssignature (x _(s)—position of the source). The source can be avolumetric source or a point-force source. The position of the source isa point in space represented by x _(s)=(x_(s), y_(s), z_(s)) where it isgenerated. The problem is considered in an unbounded domain. It isassumed that the solution satisfies a vanishing attenuation principle.The problem is solved as a first order elastic equation.

To solve Equation (1), a preconditioner to operator L is introduced asan inverse to the operator shown in Equation (13) below.

$\begin{matrix}{L_{0} = {{i\omega {M_{0}(z)}} - {P\frac{\partial}{\partial x}} - {Q\frac{\partial}{\partial y}} - {R\frac{\partial}{\partial z}}}} & (13)\end{matrix}$

In Equation (13), the matrix operator M₀ enables application of thepreconditioner in 3D rather than 2D. The matrix operator M₀ has the samestructure as the operator M shown in Equation (2) above. In addition,complex damping is introduced into the matrix operator M₀ as shown inEquations (14)-(17) below.

a ₀(z)={circumflex over (β)}â ₀(z)  (14)

b ₀(z)={circumflex over (β)}{circumflex over (b)} ₀(z)  (15)

c ₀(z)={circumflex over (β)}ĉ ₀(z)  (16)

{circumflex over (β)}=(1+iβ)  (17)

Real-valued functions, â₀(z), {circumflex over (b)}₀(z), ĉ₀(z) and ρ₀(z)are chosen. Also, β is a real constant with the value from 0 to 1. The3D operator M₀ is chosen to ensure convergence of the iterative process.Convergence depends on spectral properties of the operator L in Equation(1). To ensure convergence, the spectrum is shifted and squeezed in aproper way by implementing complex damping. Complex damping was alsoconsidered in a two-dimensional implementation. However, the operatordeveloped for the two-dimensional implementation did not improve thespectrum of the 3D operator because the 3D operator has a form that isdifferent from that of the 2D operator.

The initial operator is represented as a combination of L₀ and operatorδL as shown in Equation (18) below.

δL=−iω[M ₀(z)−M(x,y,z)]  (18)

This representation transforms Equation (1) into Equations (19) and (20)shown below.

(I−δLL ₀ ⁻¹)ũ=f  (19)

u=L ₀ ⁻¹ ũ  (20)

The computational domain is bounded to solve the system numerically. Tominimize (or avoid) any artificial reflections from the boundary of thecomputational domain, Perfectly Matched Layer (PML) is introduced in avertical direction. Alternatively, a sponge layer that is used in thelateral direction can also be introduced. A further alternative is touse first order absorbing boundary conditions (ABC). In someimplementations, a free-surface boundary condition can be chosen at thetop boundary. Introducing PML does not change the form of Equations (19)and (20). Instead, only the constant operator introduced in Equation (9)changes to R=γ(z)R, where γ(z) is a damping function. By lateraldirections we assume that the computational domain is embedded into theinfinite space or half-space filled with 1D (complex) backgroundparameters shown in equations (14)-(17), so that δL≡0 outside the targetarea. To avoid artificial reflections, the target area is encircled by atransition layer of about two wavelengths in width. In this layer, asmooth transition from the real heterogeneous medium to the backgroundmodel is performed.

The choice of depth-dependent velocity in preconditioner L₀ is chosenbased, for example, on variability in the 3D heterogeneous media. Thesubsurface is assumed to have higher velocity variability in thevertical direction than in the lateral direction, even though lateralvelocity variations exist. Therefore, the model is approximated to belaterally homogeneous. For such a model, one-dimensional (1 D)functions, â₀(z), {circumflex over (b)}₀(z), ĉ₀(z) and ρ₀(z) are chosento minimize ∥δL∥_(L) ₂ over the target area.

To avoid convergence problems related to difference in magnitude ofdifferent solution vector components (for example, u_(i)˜10⁻¹²,σ_(ij)˜10⁻⁶), the scaling shown in Equations (21) and (22) are used.

u _(i) =Ωu _(i)  (21)

σ_(ij)=Ω⁻¹σ_(ij)  (22)

In Equations (21) and (22), Ω˜10⁻³.

The iterative process is then applied to Equation (19). The solution toEquation (19) can then be multiplied by L₀ ⁻¹ to solve Equation (1). Insome implementations, the iterative process can be implemented using theKrylov-type iterative method using a matrix-free approach or othersimilar methods typically used to iteratively solve equations ingeophysical problems. For example, the bi-conjugate gradient stabilized(BiCGSTAB) method can be chosen because of moderate memory requirements.

In general, a matrix-vector multiplication process can be implemented asfast as possible. To do so, the actions of the operators δL and L₀ ⁻¹can be computed on any vector successively, instead of consideringEquation (19) as an integral equation. To compute the product L₀⁻¹u_(input), where the vector u_(input) is provided by the Krylov linearsolver, the complex dumped FD system of elasticity in 3D can be resolvedwith a depth dependent coefficient shown in Equation (23) below.

L ₀ u _(input) =u _(input)  (23)

Such resolution can be achieved by implementing two-dimensional (2D)Fast Fourier Transform (FFT) along the lateral variables followed by asolution of a number of systems of ordinary linear differentialequations with respect to depth. After the application of the Fouriertransform, several 1D systems result for each pair (k_(x), k_(y)) ofspatial frequencies. Each system is solved independently. For example, aboundary value problem for each such system is solved using afinite-different approximation of derivatives in vertical direction. Todo so, a 4^(th) order approximation staggered grid scheme can be used.The result is a set of small banded systems that can be solved quickly,for example, using an Intel Math Kernel Library. The matrix-vectorproduct procedure can be finalized using a point-by-pointmultiplication/subtraction of the given 3D functions as shown inEquation (23) below.

u _(input) =u _(input) −δL,u _(output)  (24)

The techniques described above can be parallelized via MPI using anapproach based on domain decomposition of the target area along one ofthe horizontal direction and application of the parallel version ofBiCGSTAB and FFT, as shown in FIG. 2. Doing so can enable using softwarepackages with FFT designed for modern HPC architectures to be used. EachMPI process solves its own set of small banded SLAE independently ofother MPI processes.

FIG. 1 is a flowchart of an example of a process 100 for numericalsimulation of 3D elastic waves in frequency domain. The FWI process usesthe results of the techniques described here several times at eachiteration. At 102, the 3D elastic system is represented in complianceform. At 104, corresponding medium parameters, dimensions ofcomputational domain and condition at the top boundary are set up. At106, one-dimensional (1D) background medium parameters are calculatedand preconditioner is constructed. At 108, the initial problem isrepresented as the perturbation to the preconditioner. At 110, theiterative process is started. At 112, matrix-vector multiplication isperformed at each iteration. At 114, the solution is multiplied by theinverse operator to preconditioner. The techniques described here arematrix-free. The preconditioner ensures convergence of the iterativeprocess. The techniques can be implemented using 2.5 points per minimalwavelength in lateral direction. The techniques can be extended toviscoelasticity.

Example Verification

The techniques described here were verified by comparison with an exactsolution obtained for 3D homogenous media implementing a PML boundarycondition at the top. The verification showed that the results of thetechniques described here are in very good agreement with the exactsolution.

Example Benchmarking

The techniques described here were then benchmarked with a known timedomain solver based on explicit 4^(th) order finite-difference forcomplex fragment of 3D SEG/EAGE Overthrust model with different boundaryconditions at the top (e.g., PML and free surface). The size of thecomputational domain was 200×200×187 points with the uniform grid sizeh_(x,y)=30 m in the lateral direction and h_(z)=10 m in depth. Thevertical point-force source was located in the middle of the model at 20m depth. The receivers were placed along the horizontal plane at a depthof 900 m. The wave-field was simulated for frequency of 10 Hz. Thetolerance level in the iterative process was chosen to be equal to 10⁻⁵.Some cross-sections of simulated FD u_(x)-velocity component arepresented in FIGS. 3 and 4. FIGS. 3 and 4 show plots comparing Fouriertransformed results of standard (2, 4) time domain finite-differencesimulation. The results were normalized on their maximum value. With PMLat the top boundary, the iterative process converged in 58 iterations.With free surface at the top boundary, the iterative process convergedin 80 iterations because of arising surface waves.

Example Scalability Estimation

As described above, the techniques described above can be parallelizedvia MPI. Consequently, the scalability of the parallel algorithm wastested. The results of the test are shown in FIG. 5. FIG. 5 is a plotshowing scalability of a parallel iterative solver. The scalabilitytests demonstrate the ability of the algorithm to reduce the totalcomputational time with increasing number of MPI processes and iscalculated using the formula shown in Equation (24) below.

eff _(strong)(N)=t(N)/t(N ₀)  (25)

In Equation (25), t(N) is the total computational time for N MPIprocesses, N₀ is an initial number of processes, t(N₀) is thecorresponding computational time. The homogeneous model of 800×800×200grid points with a vertical point-force source placed in the middle wasconsidered. The uniform grid step is equal to 10 m. The wave-field wasstimulated at a frequency of 5 Hz. The target area was decomposedsuccessively into a power of two number of subdomains. Each MPI processcorresponded to the whole computational node. FIG. 5 shows a scalabilitycurve in which N₀=4. The behavior is in very good agreement with theideal strong scalability curve indicating that the time spent onexchanges between MPI processes do not significantly influence theacceleration time due to parallelization.

Example Application to a Land Model

The techniques described here were applied to simulate an elasticwavefield in a typical land model having a size of about 24×21×6.5 kmwith a uniform grid step of h_(x)=h_(y)=25 m and h_(z)=10 m. Thesimulation was performed at a frequency of 5 Hz. PML boundary conditionwas used at the top surface. The vertical point-force source was placedin the middle of the target area at 20 m depth. The simulation wasperformed using 450 GB RAM implemented using 4 computational nodes with128 GB RAM each. The iterative process converged in 120 iterations. FIG.6 is a plot showing velocity distribution in a land model. FIG. 7 is aschematic diagram showing a real part of a frequency-domain horizontalvelocity component. FIG. 7 shows the real part of simulatedfrequency-domain u_(x)-velocity component (slices at Z=50 m, Y=12 km andX=9.6 km).

Implementations of the subject matter and the operations described inthis specification can be implemented in digital electronic circuitry,or in computer software, firmware, or hardware, including the structuresdisclosed in this specification and their structural equivalents, or incombinations of one or more of them. Implementations of the subjectmatter described in this specification can be implemented as one or morecomputer programs, i.e., one or more modules of computer programinstructions, encoded on computer storage medium for execution by, or tocontrol the operation of, data processing apparatus. Alternatively or inaddition, the program instructions can be encoded on anartificially-generated propagated signal, e.g., a machine-generatedelectrical, optical, or electromagnetic signal, that is generated toencode information for transmission to suitable receiver apparatus forexecution by a data processing apparatus. A computer storage medium canbe, or be included in, a computer-readable storage device, acomputer-readable storage substrate, a random or serial access memoryarray or device, or a combination of one or more of them. Moreover,while a computer storage medium is not a propagated signal, a computerstorage medium can be a source or destination of computer programinstructions encoded in an artificially-generated propagated signal. Thecomputer storage medium can also be, or be included in, one or moreseparate physical components or media (e.g., multiple CDs, disks, orother storage devices).

The operations described in this specification can be implemented asoperations performed by a data processing apparatus on data stored onone or more computer-readable storage devices or received from othersources.

The term “data processing apparatus” encompasses all kinds of apparatus,devices, and machines for processing data, including by way of example aprogrammable processor, a computer, a system on a chip, or multipleones, or combinations, of the foregoing The apparatus can includespecial purpose logic circuitry, e.g., an FPGA (field programmable gatearray) or an ASIC (application-specific integrated circuit). Theapparatus can also include, in addition to hardware, code that createsan execution environment for the computer program in question, e.g.,code that constitutes processor firmware, a protocol stack, a databasemanagement system, an operating system, a cross-platform runtimeenvironment, a virtual machine, or a combination of one or more of them.The apparatus and execution environment can realize various differentcomputing model infrastructures, such as web services, distributedcomputing and grid computing infrastructures.

A computer program (also known as a program, software, softwareapplication, script, or code) can be written in any form of programminglanguage, including compiled or interpreted languages, declarative orprocedural languages, and it can be deployed in any form, including as astand-alone program or as a module, component, subroutine, object, orother unit suitable for use in a computing environment. A computerprogram may, but need not, correspond to a file in a file system. Aprogram can be stored in a portion of a file that holds other programsor data (e.g., one or more scripts stored in a markup languagedocument), in a single file dedicated to the program in question, or inmultiple coordinated files (e.g., files that store one or more modules,sub-programs, or portions of code). A computer program can be deployedto be executed on one computer or on multiple computers that are locatedat one site or distributed across multiple sites and interconnected by acommunication network.

The processes and logic flows described in this specification can beperformed by one or more programmable processors executing one or morecomputer programs to perform actions by operating on input data andgenerating output. The processes and logic flows can also be performedby, and apparatus can also be implemented as, special purpose logiccircuitry, e.g., an FPGA (field programmable gate array) or an ASIC(application-specific integrated circuit).

Processors suitable for the execution of a computer program include, byway of example, both general and special purpose microprocessors, andany one or more processors of any kind of digital computer. Generally, aprocessor will receive instructions and data from a read-only memory ora random access memory or both. The essential elements of a computer area processor for performing actions in accordance with instructions andone or more memory devices for storing instructions and data. Generally,a computer will also include, or be operatively coupled to receive datafrom or transfer data to, or both, one or more mass storage devices forstoring data, e.g., magnetic, magneto-optical disks, or optical disks.However, a computer need not have such devices. Moreover, a computer canbe embedded in another device, e.g., a mobile telephone, a personaldigital assistant (PDA), a mobile audio or video player, a game console,a Global Positioning System (GPS) receiver, or a portable storage device(e.g., a universal serial bus (USB) flash drive), to name just a few.Devices suitable for storing computer program instructions and datainclude all forms of non-volatile memory, media and memory devices,including by way of example semiconductor memory devices, e.g., EPROM,EEPROM, and flash memory devices; magnetic disks, e.g., internal harddisks or removable disks; magneto-optical disks; and CD-ROM and DVD-ROMdisks. The processor and the memory can be supplemented by, orincorporated in, special purpose logic circuitry.

To provide for interaction with a user, implementations of the subjectmatter described in this specification can be implemented on a computerhaving a display device, e.g., a CRT (cathode ray tube) or LCD (liquidcrystal display) monitor, for displaying information to the user and akeyboard and a pointing device, e.g., a mouse or a trackball, by whichthe user can provide input to the computer. Other kinds of devices canbe used to provide for interaction with a user as well; for example,feedback provided to the user can be any form of sensory feedback, e.g.,visual feedback, auditory feedback, or tactile feedback; and input fromthe user can be received in any form, including acoustic, speech, ortactile input. In addition, a computer can interact with a user bysending documents to and receiving documents from a device that is usedby the user; for example, by sending web pages to a web browser on auser's client device in response to requests received from the webbrowser.

While this specification contains many specific implementation details,these should not be construed as limitations on the scope of anyinventions or of what may be claimed, but rather as descriptions offeatures specific to particular implementations of particularinventions. Certain features that are described in this specification inthe context of separate implementations can also be implemented incombination in a single implementation. Conversely, various featuresthat are described in the context of a single implementation can also beimplemented in multiple implementations separately or in any suitablesubcombination. Moreover, although features may be described above asacting in certain combinations and even initially claimed as such, oneor more features from a claimed combination can in some cases be excisedfrom the combination, and the claimed combination may be directed to asubcombination or variation of a subcombination.

Similarly, while operations are depicted in the drawings in a particularorder, this should not be understood as requiring that such operationsbe performed in the particular order shown or in sequential order, orthat all illustrated operations be performed, to achieve desirableresults. In certain circumstances, multitasking and parallel processingmay be advantageous. Moreover, the separation of various systemcomponents in the implementations described above should not beunderstood as requiring such separation in all implementations, and itshould be understood that the described program components and systemscan generally be integrated together in a single software product orpackaged into multiple software products.

Thus, particular implementations of the subject matter have beendescribed. Other implementations arc within the scope of the followingclaims. In some cases, the actions recited in the claims can beperformed in a different order and still achieve desirable results. Inaddition, the processes depicted in the accompanying figures do notnecessarily require the particular order shown, or sequential order, toachieve desirable results. In certain implementations, multitasking andparallel processing may be advantageous.

1-22. (canceled)
 23. A computer-implemented method comprising:representing, by data processing apparatus, a three-dimensional (3D)elastic system in compliance form, the 3D elastic system representing 3Dland media, the compliance form comprising${{Lu} = {{\left\lbrack {{i\; \omega \; {M\left( {x,y,z} \right)}} - {P\frac{\partial}{\partial x}} - {Q\frac{\partial}{\partial y}} - {R\frac{\partial}{\partial z}}} \right\rbrack u} = f}},$wherein u represents a velocity through the media, wherein M is anoperator represented by ${M = \begin{pmatrix}{\rho \; I_{3 \times 3}} & 0 \\0 & S_{6 \times 6}\end{pmatrix}},$ ρ is density, I_(3×3) is an identical operator andS_(6×6) is a matrix represented by ${S = \begin{bmatrix}a & {- b} & {- b} & 0 & 0 & 0 \\{- b} & a & {- b} & 0 & 0 & 0 \\{- b} & {- b} & a & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0\end{bmatrix}},$ a is an operator represented by${a = \frac{\lambda + \mu}{u\left( {{2\; \mu} + {3\; \lambda}} \right)}},$b is an operator represented by${b = \left( \frac{\lambda}{2\; {\mu \left( {{2\; \mu} + {3\; \lambda}} \right)}} \right)},$c is an operator represented by ${c = \left( \frac{1}{\mu} \right)},$ λand μ are Lame parameters, S_(6×6) is compliance tensor C⁻¹, P isrepresented by ${P = \begin{bmatrix}0 & A \\A^{T} & 0\end{bmatrix}},$ Q is represented by ${Q = \begin{bmatrix}0 & B \\B^{T} & 0\end{bmatrix}},$ R is represented by ${R = \begin{bmatrix}0 & D \\D^{T} & 0\end{bmatrix}},$ wherein ${A = \begin{bmatrix}1 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 1 \\0 & 0 & 0 & 0 & 1 & 0\end{bmatrix}},$ wherein ${B = \begin{bmatrix}0 & 0 & 0 & 0 & 0 & 1 \\0 & 1 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 1 & 0 & 0\end{bmatrix}},$ wherein ${D = \begin{bmatrix}0 & 0 & 0 & 0 & 1 & 0 \\0 & 0 & 0 & 1 & 0 & 0 \\0 & 0 & 1 & 0 & 0 & 0\end{bmatrix}};$ determining, by the data processing apparatus,parameters corresponding to the media, dimensions of a computationaldomain to solve the compliance form and condition at a top boundary ofthe computational domain; determining, by the data processing apparatus,one-dimensional (1D) background parameters of the media and apreconditioner solve the compliance form within the computationaldomain; representing, by the data processing apparatus, an initialproblem as a perturbation of the preconditioner; iteratively performing,by the data processing apparatus, matrix-vector multiplication to solvethe compliance form; determining, by the data processing apparatus, avelocity distribution of the velocity through the media in response toiteratively performing the matrix-vector multiplication.
 24. Thecomputer-implemented method of claim 23, wherein the 3D land mediacomprises heterogeneous media.
 25. The computer-implemented method ofclaim 23, wherein the 3D land media comprises isotropic media.
 26. Thecomputer-implemented method of claim 23, wherein the 3D land mediacomprises a rock formation.
 27. The computer-implemented method of claim23, wherein the 3D land media comprises a hydrocarbon reservoir.
 28. Thecomputer-implemented method of claim 23, wherein determining thepreconditioner comprises constructing the preconditioner as an inverseto the operator shown in the following equation:${L_{0} = \; {{i\; \omega \; {M_{0}(z)}} - {P\frac{\partial}{\partial x}} - {Q\frac{\partial}{\partial y}} - {R\frac{\partial}{\partial z}}}},$wherein the operator M₀ has the same structure as the operator M. 29.The computer-implemented method of claim 28, wherein a complex dampingfunction is introduced into the operator M₀.
 30. Thecomputer-implemented method of claim 29, wherein the complex dampingfunction is represented by the following equations: a₀(z)={circumflexover (β)}â₀(z), b₀(z)={circumflex over (β)}{circumflex over (b)}₀(z),c₀(z)={circumflex over (β)}ĉ₀(z), {circumflex over (β)}=(1+iβ).
 31. Thecomputer-implemented method of claim 23, wherein representing theinitial problem as the perturbation to the preconditioner comprisesrepresenting the initial operator L as a combination of L₀ and δL,wherein representing the initial operator L as a combination of L₀ andδL comprises representing the initial operator according to thefollowing equation: δL=−iω[M₀(z)−M(x,y,z)].
 32. The computer-implementedmethod of claim 23, wherein the matrix-vector multiplication isperformed by the application of a parallel Fast Fourier Transform (FFT),a banded matrix solver, and an inverse FFT.